Indicador Social en la niñez Colombiana usando Machine Learning
Beatriz Elena Jaramillo Gallego
20/12/2020
1. Data Set (ENCV-2018)
Está parte del proyecto fue desarrollada en Stata.
La Encuesta de Calidad de Vida (ECV) es una investigación que el DANE (Departamento Administrativo Nacional de Estadística) realiza con el objeto de recoger información sobre diferentes aspectos y dimensiones del bienestar y las condiciones de vida de los hogares, incluyendo temas como: el acceso a bienes y servicios públicos, privados o comunales, salud, educación, atención integral de niños y niñas menores de 5 años, entre otros. (https://www.dane.gov.co/index.php/estadisticas-por-tema/salud/calidad-de-vida-ecv/encuesta-nacional-de-calidad-de-vida-ecv-2018)
2. Análisis exploratorio
Se hizo un análisis exploratorio descriptivo de la base de datos de las variables que se han seleccionado para continuar con el análisis que se ha propuesto para este trabajo de final de Master.
Tabla resumen Descripción general amplia de un dataframe.
# La función skim() es una alternativa a summary(), la cual nos proporciona rápidamente una descripción general amplia de un dataframe.
skim(encv_2018_nens) | Name | encv_2018_nens |
| Number of rows | 26957 |
| Number of columns | 27 |
| _______________________ | |
| Column type frequency: | |
| factor | 20 |
| numeric | 7 |
| ________________________ | |
| Group variables | None |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| Sexo | 0 | 1 | FALSE | 2 | Hom: 13843, Muj: 13114 |
| vacunas | 0 | 1 | FALSE | 2 | Si: 25478, No: 1479 |
| control_crec | 0 | 1 | FALSE | 2 | Si: 23698, No: 3259 |
| icbf | 0 | 1 | FALSE | 2 | No: 25423, Si: 1534 |
| ingresos_flia | 0 | 1 | FALSE | 3 | No : 13541, Sol: 12420, Cub: 996 |
| seg_alim | 0 | 1 | FALSE | 2 | No: 25160, Si: 1797 |
| REGION | 0 | 1 | FALSE | 9 | Car: 8139, Ori: 5459, Cen: 4453, Ori: 4146 |
| Zona | 0 | 1 | FALSE | 2 | Urb: 13757, Rur: 13200 |
| Registro_civil | 0 | 1 | FALSE | 2 | Si: 25772, No: 1185 |
| hacinnomiti | 0 | 1 | FALSE | 2 | Sin: 24727, Con: 2230 |
| quintil | 0 | 1 | FALSE | 5 | Riq: 5570, Ric: 5537, Muy: 5406, Pob: 5377 |
| sgss | 0 | 1 | FALSE | 2 | Si: 24766, No: 2191 |
| Estado_de_Salud | 0 | 1 | FALSE | 2 | Muy: 25610, Reg: 1347 |
| No_Actividad_Estimulacion | 0 | 1 | FALSE | 2 | No: 17267, Si: 9690 |
| Actividad_Deportiva | 0 | 1 | FALSE | 2 | No: 26050, Si: 907 |
| Victima_Hecho_Viol | 0 | 1 | FALSE | 2 | No: 24405, Si: 2552 |
| Afect_Fenome_Natural | 0 | 1 | FALSE | 2 | No: 22166, Si: 4791 |
| cuidado_menor | 0 | 1 | FALSE | 2 | No: 26887, Si: 70 |
| hogar_exclusion | 0 | 1 | FALSE | 2 | No: 15015, Si: 11942 |
| DEPTO_D | 0 | 1 | FALSE | 33 | Gua: 1459, Ces: 1186, Atl: 1185, Bol: 1142 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| DIRECTORIO | 0 | 1 | 7058700.77 | 33185.39 | 7000011.00 | 7029096.00 | 7060112.00 | 7086896.00 | 7119632.00 | <U+2587><U+2586><U+2587><U+2587><U+2586> |
| SECUENCIA_ENCUESTA | 0 | 1 | 4.36 | 1.64 | 2.00 | 3.00 | 4.00 | 5.00 | 20.00 | <U+2587><U+2582><U+2581><U+2581><U+2581> |
| SECUENCIA_P | 0 | 1 | 1.01 | 0.13 | 1.00 | 1.00 | 1.00 | 1.00 | 4.00 | <U+2587><U+2581><U+2581><U+2581><U+2581> |
| FEX_C | 0 | 1 | 186.62 | 316.98 | 2.43 | 55.35 | 107.53 | 208.77 | 3356.88 | <U+2587><U+2581><U+2581><U+2581><U+2581> |
| Edad | 0 | 1 | 2.59 | 1.70 | 0.00 | 1.00 | 3.00 | 4.00 | 5.00 | <U+2587><U+2585><U+2585><U+2585><U+2585> |
| I_HOGAR | 0 | 1 | 1592578.88 | 2909244.14 | 0.00 | 584166.67 | 1040000.00 | 1799833.33 | 313045000.00 | <U+2587><U+2581><U+2581><U+2581><U+2581> |
| PERCAPITA | 0 | 1 | 361400.79 | 587206.45 | 0.00 | 125000.00 | 230000.00 | 395625.00 | 34782777.78 | <U+2587><U+2581><U+2581><U+2581><U+2581> |
Distribución niños por departamentos
encv_2018_nens %>% tab(DEPTO_D)##
## DEPTO_D ¦ Freq. Percent Cum.
## -------------------------+----------------------------
## Amazonas ¦ 758 2.81 2.81
## Antioquia ¦ 784 2.91 5.72
## Arauca ¦ 620 2.30 8.02
## Atlantico ¦ 1185 4.40 12.42
## Bogota ¦ 603 2.24 14.65
## Bolivar ¦ 1142 4.24 18.89
## Boyaca ¦ 578 2.14 21.03
## Caldas ¦ 620 2.30 23.33
## Caqueta ¦ 936 3.47 26.81
## Casanare ¦ 771 2.86 29.67
## Cauca ¦ 892 3.31 32.97
## Cesar ¦ 1186 4.40 37.37
## Choco ¦ 1055 3.91 41.29
## Cordoba ¦ 987 3.66 44.95
## Cundinamarca ¦ 803 2.98 47.93
## Guainia ¦ 549 2.04 49.96
## Guajira ¦ 1459 5.41 55.38
## Guaviare ¦ 636 2.36 57.74
## Huila ¦ 913 3.39 61.12
## Magdalena ¦ 1131 4.20 65.32
## Meta ¦ 910 3.38 68.69
## Narino ¦ 788 2.92 71.62
## Norte de Santander ¦ 829 3.08 74.69
## Putumayo ¦ 689 2.56 77.25
## Quindio ¦ 613 2.27 79.52
## Risaralda ¦ 605 2.24 81.77
## San Andres y Providencia ¦ 235 0.87 82.64
## Santander ¦ 726 2.69 85.33
## Sucre ¦ 1049 3.89 89.22
## Tolima ¦ 766 2.84 92.07
## Valle del Cauca ¦ 703 2.61 94.67
## Vaupes ¦ 731 2.71 97.38
## Vichada ¦ 705 2.62 100.00Distribución niños por Edad y Sexo
encv_2018_nens %>% tab(Edad,Sexo)##
## Edad ¦ Sexo ¦ Freq. Percent Cum.
## -----+--------+----------------------------
## 0 ¦ Hombre ¦ 2065 7.66 7.66
## 0 ¦ Mujer ¦ 2010 7.46 15.12
## -----+--------+----------------------------
## 1 ¦ Hombre ¦ 2296 8.52 23.63
## 1 ¦ Mujer ¦ 2101 7.79 31.43
## -----+--------+----------------------------
## 2 ¦ Hombre ¦ 2336 8.67 40.09
## 2 ¦ Mujer ¦ 2115 7.85 47.94
## -----+--------+----------------------------
## 3 ¦ Hombre ¦ 2316 8.59 56.53
## 3 ¦ Mujer ¦ 2206 8.18 64.71
## -----+--------+----------------------------
## 4 ¦ Hombre ¦ 2403 8.91 73.63
## 4 ¦ Mujer ¦ 2329 8.64 82.27
## -----+--------+----------------------------
## 5 ¦ Hombre ¦ 2427 9.00 91.27
## 5 ¦ Mujer ¦ 2353 8.73 100.00Número de datos ausentes por variable
encv_2018_nens %>% map_dbl(.f = function(x){sum(is.na(x))})## DIRECTORIO SECUENCIA_ENCUESTA SECUENCIA_P
## 0 0 0
## FEX_C Sexo Edad
## 0 0 0
## vacunas control_crec icbf
## 0 0 0
## I_HOGAR PERCAPITA ingresos_flia
## 0 0 0
## seg_alim REGION Zona
## 0 0 0
## Registro_civil hacinnomiti quintil
## 0 0 0
## sgss Estado_de_Salud No_Actividad_Estimulacion
## 0 0 0
## Actividad_Deportiva Victima_Hecho_Viol Afect_Fenome_Natural
## 0 0 0
## cuidado_menor hogar_exclusion DEPTO_D
## 0 0 0Las variables seleccionadas de la base de datos no tienen valores ausentes.
3. Preprocesamiento
Se realizarán las transformaciones necesarias en la búsqueda de encontrar buenos resultados en el algoritmo de machine learning. El modelo será aprendido con los datos de entrenamiento (80%) y luego se aplicará el modelo a los datos de test (20%), luego se evaluarán los modelos que se generen y el error incurrido.
Imputación de valores ausentes
nrow(na.omit(encv_2018_ml))## [1] 26957No tenemos valores ausentes en las variables seleccionadas para nuestro estudio.
Binarización de las variables cualitativas
Se crean variables dummy con cada uno de los niveles de las variables cualitativas de las variables seleccionadas anteriormente.
head(encv_2018_ml)| FEX_C | Sexo | Edad | vacunas | control_crec | icbf | I_HOGAR | ingresos_flia | seg_alim | Zona | Registro_civil | hacinnomiti | quintil | sgss | Estado_de_Salud | No_Actividad_Estimulacion | Actividad_Deportiva | Victima_Hecho_Viol | Afect_Fenome_Natural | cuidado_menor | hogar_exclusion | DEPTO_D |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 468.5849 | 2 | 1 | 1 | 1 | 2 | 233833.3 | 1 | 2 | 2 | 1 | 2 | 2 | 1 | 1 | 2 | 2 | 2 | 2 | 2 | Si | 2 |
| 353.7918 | 2 | 2 | 1 | 1 | 2 | 976666.7 | 2 | 2 | 2 | 1 | 2 | 3 | 1 | 1 | 2 | 2 | 2 | 2 | 2 | No | 2 |
| 399.9132 | 2 | 0 | 1 | 1 | 2 | 759166.7 | 1 | 2 | 2 | 1 | 2 | 1 | 1 | 1 | 1 | 2 | 2 | 2 | 2 | Si | 2 |
| 399.9132 | 1 | 5 | 1 | 1 | 2 | 475000.0 | 1 | 1 | 2 | 1 | 2 | 3 | 1 | 1 | 1 | 2 | 2 | 2 | 2 | No | 2 |
| 921.8979 | 1 | 2 | 1 | 1 | 2 | 993308.7 | 2 | 2 | 1 | 1 | 2 | 5 | 1 | 1 | 1 | 2 | 2 | 2 | 2 | No | 2 |
| 978.2812 | 2 | 4 | 1 | 1 | 2 | 2224666.7 | 1 | 2 | 1 | 1 | 2 | 3 | 1 | 1 | 2 | 1 | 1 | 2 | 2 | No | 2 |
Variables con varianza próxima a cero
DT::datatable(binari,
extensions = list('Scroller'=TRUE ,'FixedColumns'=NULL),
options = list(autowidth=T,pageLength=20,fixedHeader=T,scrollY=200,scrollX=200,searching = T,
fixedColumns=list(leftColumns=1) )) %>% formatRound('freqRatio', 3) %>% formatPercentage('percentUnique', 2)Entre los predictores incluidos en el modelo, se detecta que los predictores: Registro_civil, Actividad_Deportiva, cuidado_menor y Estado_de_Salud cuentan con varianza cero o próxima a cero.
Estandarización y escalado
La escala cuando los predictores son numéricos, así como la magnitud de su varianza pueden influir en el modelo. He usado la estrategia de centrado en la variable Ingresos del hogar, para evitar que está tenga una influencia en el modelo.
encv_2018_ml$I_HOGAR <- scale(encv_2018_ml$I_HOGAR, center=TRUE, scale=FALSE)
encv_2018_ml$I_HOGAR <- as.numeric(encv_2018_ml$I_HOGAR)
head(encv_2018_ml$I_HOGAR)## [1] -1358745.5 -615912.2 -833412.2 -1117578.9 -599270.2 632087.8Distribución de variables respuesta
Cuando se crea un modelo, es muy importante estudiar la distribución de la variable respuesta, ya que, a fin de cuentas, es lo que nos interesa predecir.
ggplot(data = encv_2018_ml, aes(x = hogar_exclusion, y = ..count.., fill = hogar_exclusion)) +
geom_bar() +
scale_fill_manual(values = c("gray50", "orangered2")) +
labs(title = "Hogares en Exclusión") +
theme_bw() +
theme(legend.position = "bottom")
# Tabla de frecuencias
table(encv_2018_ml$hogar_exclusion)| Si | No |
|---|---|
| 11942 | 15015 |
prop.table(table(encv_2018_ml$hogar_exclusion)) %>% round(digits = 2)| Si | No |
|---|---|
| 0.44 | 0.56 |
Porcentaje de aciertos si se predicen los hogares con exclusión social.
El porcentaje mínimo que hay que intentar superar clasificar a los hogares con exclusion social con los modelos predictivos es del 55.7%.
n_observaciones <- nrow(encv_2018_ml)
predicciones <- rep(x = "No", n_observaciones)
mean(predicciones == encv_2018_ml$hogar_exclusion) * 100## [1] 55.69982Distribución de la variable respuesta vs nivel de riqueza
ggplot(data = encv_2018_ml, aes(x = quintil, y = ..count.., fill = hogar_exclusion)) +
geom_bar() +
scale_fill_manual(values = c("gray50", "orangered2")) +
labs(title = "Hogares en Exclusión VS Nivel de Riqueza") +
theme_bw() +
theme(legend.position = "bottom")
Distribución de la variable respuesta vs Hacinamiento
ggplot(data = encv_2018_ml, aes(x = hacinnomiti, y = ..count.., fill = hogar_exclusion)) +
geom_bar() +
scale_fill_manual(values = c("gray50", "orangered2")) +
labs(title = "Hogares en Exclusión VS Hacinamiento") +
theme_bw() +
theme(legend.position = "bottom")
encv_2018_nens %>% tab(hacinnomiti,hogar_exclusion)##
## hacinnomiti ¦ hogar_exclusion ¦ Freq. Percent Cum.
## -----------------+-----------------+----------------------------
## Con hacinamiento ¦ Si ¦ 1854 6.88 6.88
## Con hacinamiento ¦ No ¦ 376 1.39 8.27
## -----------------+-----------------+----------------------------
## Sin Hacinamiento ¦ Si ¦ 10088 37.42 45.69
## Sin Hacinamiento ¦ No ¦ 14639 54.31 100.00División de los datos en entrenamiento y test
La base de datos de niños se ha divido en 80% para “Training” y 20% para “Test”, estos porcentajes suele dar buenos resultados. Esta repartición se ha realizado de forma aleatoria-estratificada.
Se genera una base de datos en donde se ha excluido las variables DIRECTORIO, SECUENCIA_ENCUESTA, SECUENCIA_P, FEX_C, PERCAPITA, REGION, SEXO, CUIDADO_MENOR,-ESTADO_DE_SALUD, VICTIMA_HECHO_VIOL, VACUNAS y EDAD que no son relevantes en nuestro análisis posterior y que solo se han usado para tener un orden en el momento de generar la base de datos de niños.
Se verifica a continuación que la distribución de la variable respuesta sea similar tanto en el conjunto de entrenamiento y como en el conjunto de test.
Training
round(prop.table(table(datos_train$hogar_exclusion))*100,3)| Si | No |
|---|---|
| 44.301 | 55.699 |
Test
round(prop.table(table(datos_test$hogar_exclusion))*100,3)| Si | No |
|---|---|
| 44.296 | 55.704 |
Este reparto estratificado de la base de datos asegura que el conjunto de entrenamiento y el conjunto de test sean similares en cuanto a la variable respuesta.
glimpse(datos_train)## Rows: 21,566
## Columns: 22
## $ FEX_C <dbl> 468.5849, 399.9132, 399.9132, 921.8979, 9...
## $ Sexo <int> 2, 2, 1, 1, 2, 1, 2, 2, 2, 1, 1, 1, 1, 1,...
## $ Edad <int> 1, 0, 5, 2, 4, 1, 4, 5, 5, 4, 4, 1, 2, 5,...
## $ vacunas <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ control_crec <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1,...
## $ icbf <int> 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2,...
## $ I_HOGAR <dbl> -1358745.54, -833412.21, -1117578.88, -59...
## $ ingresos_flia <int> 1, 1, 1, 2, 1, 2, 1, 2, 1, 2, 1, 1, 2, 2,...
## $ seg_alim <int> 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2,...
## $ Zona <int> 2, 2, 2, 1, 1, 2, 2, 1, 2, 1, 2, 2, 2, 1,...
## $ Registro_civil <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1,...
## $ hacinnomiti <int> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2,...
## $ quintil <int> 2, 1, 3, 5, 3, 4, 3, 5, 3, 3, 2, 2, 4, 4,...
## $ sgss <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1,...
## $ Estado_de_Salud <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ No_Actividad_Estimulacion <int> 2, 1, 1, 1, 2, 2, 2, 1, 1, 2, 1, 2, 1, 1,...
## $ Actividad_Deportiva <int> 2, 2, 2, 2, 1, 2, 2, 2, 2, 1, 2, 2, 2, 2,...
## $ Victima_Hecho_Viol <int> 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2,...
## $ Afect_Fenome_Natural <int> 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2,...
## $ cuidado_menor <int> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,...
## $ hogar_exclusion <fct> Si, Si, No, No, No, No, No, No, No, No, S...
## $ DEPTO_D <int> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,...glimpse(datos_test)## Rows: 5,391
## Columns: 22
## $ FEX_C <dbl> 353.7918, 981.3012, 596.4729, 303.5955, 3...
## $ Sexo <int> 2, 2, 1, 2, 1, 1, 2, 1, 1, 2, 1, 2, 1, 1,...
## $ Edad <int> 2, 3, 1, 5, 2, 1, 3, 0, 0, 4, 5, 2, 2, 3,...
## $ vacunas <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ control_crec <int> 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1,...
## $ icbf <int> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,...
## $ I_HOGAR <dbl> -615912.2, -811578.9, 11455087.8, -118757...
## $ ingresos_flia <int> 2, 1, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2,...
## $ seg_alim <int> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,...
## $ Zona <int> 2, 1, 2, 2, 2, 1, 1, 1, 2, 1, 1, 1, 1, 2,...
## $ Registro_civil <int> 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ hacinnomiti <int> 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2,...
## $ quintil <int> 3, 5, 1, 1, 3, 5, 4, 2, 5, 3, 4, 5, 5, 1,...
## $ sgss <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ Estado_de_Salud <int> 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1,...
## $ No_Actividad_Estimulacion <int> 2, 1, 2, 1, 2, 2, 2, 2, 2, 2, 1, 2, 2, 1,...
## $ Actividad_Deportiva <int> 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2,...
## $ Victima_Hecho_Viol <int> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2,...
## $ Afect_Fenome_Natural <int> 2, 2, 2, 1, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2,...
## $ cuidado_menor <int> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,...
## $ hogar_exclusion <fct> No, No, Si, Si, No, No, No, Si, No, No, N...
## $ DEPTO_D <int> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,...Tras el Preprocesado de los datos, se han generado un total de 21 variables (20 predictores y la variable respuesta).
Selección de Predictores
Se cuenta con 21 variables que se incluiran como predictores ya que están realmente relacionadas con la variable respuesta, ya que son las que contienen información útil para el modelo, ya que incluir un exceso de variables suele conllevar una reducción de la capacidad predictiva del modelo cuando se expone a nuevos datos (overfitting).
Métodos de filtrado
Los métodos basados en filtrado evalúan la relevancia de los predictores fuera del modelo para, posteriormente, incluir únicamente aquellos predictores con un p-value inferior a un determinado límite o los n mejores.
rf_sbf##
## Selection By Filter
##
## Outer resampling method: Cross-Validated (10 fold, repeated 1 times)
##
## Resampling performance:
##
## Accuracy Kappa AccuracySD KappaSD
## 0.9998 0.9996 0.000391 0.0007925
##
## Using the training set, 18 variables were selected:
## Edad, vacunas, control_crec, icbf, I_HOGAR...
##
## During resampling, the top 5 selected variables (out of a possible 18):
## Actividad_Deportiva (100%), Afect_Fenome_Natural (100%), control_crec (100%), Edad (100%), Estado_de_Salud (100%)
##
## On average, 17.8 variables were selected (min = 17, max = 18)rf_sbf$optVariables## [1] "Edad" "vacunas"
## [3] "control_crec" "icbf"
## [5] "I_HOGAR" "ingresos_flia"
## [7] "seg_alim" "Zona"
## [9] "Registro_civil" "hacinnomiti"
## [11] "quintil" "sgss"
## [13] "Estado_de_Salud" "No_Actividad_Estimulacion"
## [15] "Actividad_Deportiva" "Victima_Hecho_Viol"
## [17] "Afect_Fenome_Natural" "DEPTO_D"La selección del método de filtrado identifica como mejor modelo el formado por 18 predictores.
Importancia de las Variables
Otra estrategia ampliamente extendida para estudiar la importancia de variables es el empleo de Random Forest.
ggarrange(p1, p2)
Ambos análisis apuntan a que las variables quintil, icbf, Registro Civil, hacinamiento, I_HOGRA tienen una influencia alta sobre las probabilidades de Exclusión Social.
4. Creación de un Modelo Predictivo
Los datos han sido Preprocesados y se han seleccionado los predictores, lo siguiente es emplear un algoritmo de machine learning que nos permita crear un modelo capaz de representar los patrones en los datos de entrenamiento y testearlos en nuevas observaciones. Por lo general, las etapas seguidas para obtener un buen modelo son:
- Entrenar el modelo
- Predicción y evaluación del modelo
- Mejora del modelo
Entrenamiento del modelo
Entre los argumentos de la función train() que se plantean usar se destacan:
* formula: la fórmula del modelo que se quiere crear.
* x, y: se pueden pasar por separado los valores de los predictores y de la variable respuesta.
* method: el nombre del algoritmo que se desea emplear.
* metric: las métricas empleadas para evaluar la capacidad predictiva del modelo.
* trControl: especificaciones adicionales sobre la forma de llevar a cabo el entrenamiento del modelo.
* …: argumentos propios del algoritmo empleado.
Ajustaré un modelo basado en una máquina vector soporte lineal que prediga la exclusión social en Colombia por departamentos en función de los predictores disponibles, pero tambien se excluirán las variables con poca importancia o influencia en la creación del modelo predictivo.
# Predictores: ## "control_crec", "icbf","I_HOGAR", "ingresos_flia", "seg_alim", "Registro_civil", "hacinnomiti", "quintil", "sgss" , "Actividad_Deportiva","Afect_Fenome_Natural","Zona","DEPTO_D"
datos_train <- select(datos_train, -Sexo, -cuidado_menor,-Estado_de_Salud,-Victima_Hecho_Viol,-vacunas,-Edad)
modelo_svmlineal <- train(hogar_exclusion ~ .,
method = "svmLinear",
data = datos_train %>% select(-FEX_C))modelo_svmlineal$finalModel## Support Vector Machine object of class "ksvm"
##
## SV type: C-svc (classification)
## parameter : cost C = 1
##
## Linear (vanilla) kernel function.
##
## Number of Support Vectors : 1042
##
## Objective Function Value : -719.1111
## Training error : 0.012566Evaluación del modelo mediante resampling
A continuación se ajusta de nuevo una máquina vector soporte lineal, esta vez con validación cruzada repetida para estimar su error.
modelo_svmlineal## Support Vector Machines with Linear Kernel
##
## 21566 samples
## 14 predictor
## 2 classes: 'Si', 'No'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 1 times)
## Summary of sample sizes: 19409, 19409, 19408, 19410, 19410, 19410, ...
## Resampling results:
##
## Accuracy Kappa
## 0.9874342 0.97453
##
## Tuning parameter 'C' was held constant at a value of 1DT::datatable(modelo_svmlineal$resample,
extensions = list('Scroller'=TRUE ,'FixedColumns'=NULL),
options = list(autowidth=T,pageLength=20,fixedHeader=T,scrollY=200,scrollX=200,searching = T,
fixedColumns=list(leftColumns=1) )) %>% formatPercentage(1:2, 3)final_plot
El accuracy promedio estimado usando un modelo de Máquina Vector Soporte Lineal, se ha conseguido un accuracy promedio de validación del 98.74%, esto significa que el modelo predice correctamente la exclusion social por departamentos para Colombia un 98.74% de las veces.
Optimización de hiperparámetros
modelo_svmlineal## Support Vector Machines with Linear Kernel
##
## 21566 samples
## 14 predictor
## 2 classes: 'Si', 'No'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 3 times)
## Summary of sample sizes: 19409, 19409, 19408, 19410, 19410, 19410, ...
## Resampling results across tuning parameters:
##
## C Accuracy Kappa
## 1e-03 0.9840952 0.9677457
## 1e-02 0.9874338 0.9745306
## 1e-01 0.9874338 0.9745306
## 5e-01 0.9874338 0.9745306
## 1e+00 0.9874338 0.9745306
## 1e+01 0.9874338 0.9745306
##
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was C = 0.01.Como hemos empleado validación cruzada con 10 particiones y 1 repeticiones, caret ha ajustado el modelo svmlineal 10 x 3 = 30 veces por cada valor de C.
DT::datatable(modelo_svmlineal$resample,
extensions = list('Scroller'=TRUE ,'FixedColumns'=NULL),
options = list(autowidth=T,pageLength=20,fixedHeader=T,scrollY=200,scrollX=200,searching = T,
fixedColumns=list(leftColumns=1) )) %>% formatPercentage(1:2, 3)ggplot(data = modelo_svmlineal$resample,
aes(x = as.factor(C), y = Accuracy, color = as.factor(C))) +
geom_boxplot(outlier.shape = NA, alpha = 0.6) +
geom_jitter(width = 0.2, alpha = 0.6) +
# Línea horizontal en el accuracy basal
geom_hline(yintercept = 0.98, linetype = "dashed") +
labs(x = "C") +
theme_bw() + theme(legend.position = "none")
ggplot(modelo_svmlineal, highlight = TRUE) +
labs(title = "Evolución del accuracy del modelo_svmlineal en función de C") +
theme_bw()
5. Predicción
Se ajustara el modelo usando la función predict() del paquete caret, con el que podremos predecir nuevos datos utilizando el objeto devuelto por train().
predicciones_raw_svmlin <- predict(modelo_svmlineal,
newdata = datos_test %>% select(-FEX_C),
type = "raw")head(predicciones_raw_svmlin)## [1] No No Si Si No No
## Levels: Si NoEl algoritmo svmLinear no calcula probabilidades, para obtenerlas se reajustara el modelo indicando en el control_train classProbs = TRUE.
predicciones_prob_svml <- predict(modelo_svmlineal,
newdata = datos_test %>% select(-FEX_C),
type = "prob")
predicciones_prob_svml %>% head()| Si | No |
|---|---|
| 0.0564246 | 0.9435754 |
| 0.0000006 | 0.9999994 |
| 0.9998236 | 0.0001764 |
| 1.0000000 | 0.0000000 |
| 0.0564246 | 0.9435754 |
| 0.0000006 | 0.9999994 |
testY = datos_test %>% select(-FEX_C)
predicciones <- extractPrediction(
models = list(svm = modelo_svmlineal),
testX = datos_test %>% select(-hogar_exclusion,-FEX_C),
testY = testY$hogar_exclusion )## Warning in method$predict(modelFit = modelFit, newdata = newdata, submodels =
## param): kernlab class prediction calculations failed; returning NAspredicciones %>% head()| obs | pred | model | dataType | object |
|---|---|---|---|---|
| Si | Si | svmLinear | Training | svm |
| Si | Si | svmLinear | Training | svm |
| No | No | svmLinear | Training | svm |
| No | No | svmLinear | Training | svm |
| No | No | svmLinear | Training | svm |
| No | No | svmLinear | Training | svm |
Error de test
confusionMatrix(data = predicciones_raw_svmlin,
reference = datos_test$hogar_exclusion,
positive = "Si")## Confusion Matrix and Statistics
##
## Reference
## Prediction Si No
## Si 2349 18
## No 39 2985
##
## Accuracy : 0.9894
## 95% CI : (0.9863, 0.992)
## No Information Rate : 0.557
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.9786
##
## Mcnemar's Test P-Value : 0.008071
##
## Sensitivity : 0.9837
## Specificity : 0.9940
## Pos Pred Value : 0.9924
## Neg Pred Value : 0.9871
## Prevalence : 0.4430
## Detection Rate : 0.4357
## Detection Prevalence : 0.4391
## Balanced Accuracy : 0.9888
##
## 'Positive' Class : Si
## # Error de test
error_test_svmlin <- mean(predicciones_raw_svmlin != datos_test$hogar_exclusion)
paste("El error del test:", round(error_test_svmlin*100, 2), "%")## [1] "El error del test: 1.06 %"6. Modelos
A continuación se entrenan diferentes modelos de machine learning con el objetivo de compararlos e identificar cual de estos arroja el mejor resultado obtiene prediciendo la exclusión social en Colombia por departamentos.
K-Nearest Neighbor (kNN)
K-Nearest Neighbor es uno de los algoritmos de machine learning más simples, este algoritmo intentará identificar observaciones en el conjunto de entrenamiento que se asemejen a la observación de test y asignarle como valor predicho la clase predominante entre dichas observaciones.
Ajuste, optimización y validación del modelo
Se ajusta el K-Nearest Neighbor, esta vez con validación cruzada repetida para estimar su error
__*k:__ número de observaciones vecinas empleadas.
predicciones_knn <- predict(modelo_knn,
newdata = datos_train %>% select(-FEX_C))modelo_knn## k-Nearest Neighbors
##
## 21566 samples
## 14 predictor
## 2 classes: 'Si', 'No'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 5 times)
## Summary of sample sizes: 19409, 19409, 19408, 19410, 19410, 19410, ...
## Resampling results across tuning parameters:
##
## k Accuracy Kappa
## 1 0.7753412 0.5451179
## 2 0.7257998 0.4446957
## 5 0.7218215 0.4351415
## 10 0.7088847 0.4075520
## 15 0.7102666 0.4101846
## 20 0.7024766 0.3941458
## 30 0.6958175 0.3805619
## 50 0.6930539 0.3747924
##
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was k = 1.Con un modelo kNN con k=1 se consigue un accuracy de validación promedio del 77.53%
# Gráfica del modelo
# ==============================================================================
ggplot(modelo_knn, highlight = TRUE) +
scale_x_continuous(breaks = hiperparametros$k) +
labs(title = "Evolución del accuracy del modelo KNN", x = "K") +
theme_bw()
summary(modelo_knn$resample$Accuracy)| Min. | 1st Qu. | Median | Mean | 3rd Qu. | Max. |
|---|---|---|---|---|---|
| 0.6675939 | 0.6993275 | 0.7106422 | 0.7166827 | 0.7244898 | 0.788961 |
El accuracy promedio estimado usando un modelo de K-Nearest Neighbor, se ha conseguido un accuracy promedio de validación del 71.67%, esto significa que el modelo predice correctamente la exclusion social por departamentos para Colombia un 71.67% de las veces.
Error de test
confusionMatrix(data = predicciones_raw_knn,
reference = datos_test$hogar_exclusion,
positive = "Si")## Confusion Matrix and Statistics
##
## Reference
## Prediction Si No
## Si 1832 603
## No 556 2400
##
## Accuracy : 0.785
## 95% CI : (0.7738, 0.7959)
## No Information Rate : 0.557
## P-Value [Acc > NIR] : <2e-16
##
## Kappa : 0.5652
##
## Mcnemar's Test P-Value : 0.1766
##
## Sensitivity : 0.7672
## Specificity : 0.7992
## Pos Pred Value : 0.7524
## Neg Pred Value : 0.8119
## Prevalence : 0.4430
## Detection Rate : 0.3398
## Detection Prevalence : 0.4517
## Balanced Accuracy : 0.7832
##
## 'Positive' Class : Si
## # Error de test
error_test_KNN <- mean(predicciones_raw_knn != datos_test$hogar_exclusion)
paste("El error del test:", round(error_test_KNN*100, 2), "%")## [1] "El error del test: 21.5 %"Regresión logística
Ajuste, optimización y validación del modelo
Este algoritmo no tiene ningún hiperparámetro pero, para que efectúe una regresión logística, hay que indicar family = “binomial”.
predicciones_logistic <- predict(modelo_logistic,
newdata = datos_train %>% select(-FEX_C))Una validación cruzada con 10 particiones y 5 repeticiones implica ajustar y evaluar el modelo 10 x 5 = 50 veces, cada vez con una partición distinta, y un ajuste con los datos de entrenamiento para crear el modelo final.
modelo_logistic## Generalized Linear Model
##
## 21566 samples
## 14 predictor
## 2 classes: 'Si', 'No'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 5 times)
## Summary of sample sizes: 19409, 19409, 19408, 19410, 19410, 19410, ...
## Resampling results:
##
## Accuracy Kappa
## 0.9871001 0.9738518summary(modelo_logistic$finalModel)##
## Call:
## NULL
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -5.3688 -0.0115 0.0006 0.0164 4.4829
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.357e+00 2.057e+02 -0.026 0.979218
## control_crec -1.380e-01 2.049e-01 -0.674 0.500609
## icbf 1.287e+01 3.243e-01 39.676 < 2e-16 ***
## I_HOGAR -1.127e-07 3.062e-08 -3.680 0.000234 ***
## ingresos_flia -5.291e-02 1.000e-01 -0.529 0.596785
## seg_alim 8.114e-01 2.003e-01 4.052 5.08e-05 ***
## Zona -4.935e-02 1.049e-01 -0.471 0.637918
## Registro_civil -3.158e+01 2.057e+02 -0.154 0.877954
## hacinnomiti 1.057e+00 2.061e-01 5.127 2.95e-07 ***
## quintil 6.266e+00 1.219e-01 51.406 < 2e-16 ***
## sgss 1.706e-02 2.405e-01 0.071 0.943462
## No_Actividad_Estimulacion -1.052e+00 1.214e-01 -8.663 < 2e-16 ***
## Actividad_Deportiva -3.198e+00 2.123e-01 -15.061 < 2e-16 ***
## Afect_Fenome_Natural 8.927e-02 1.267e-01 0.705 0.481058
## DEPTO_D -4.156e-03 5.504e-03 -0.755 0.450119
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 29616.1 on 21565 degrees of freedom
## Residual deviance: 3145.4 on 21551 degrees of freedom
## AIC: 3175.4
##
## Number of Fisher Scoring iterations: 18DT::datatable(modelo_logistic$resample,
extensions = list('Scroller'=TRUE ,'FixedColumns'=NULL),
options = list(autowidth=T,pageLength=20,fixedHeader=T,scrollY=200,scrollX=200,searching = T,
fixedColumns=list(leftColumns=1) )) %>% formatPercentage(1:2, 3)summary(modelo_logistic$resample$Accuracy)| Min. | 1st Qu. | Median | Mean | 3rd Qu. | Max. |
|---|---|---|---|---|---|
| 0.9809833 | 0.9860853 | 0.9865616 | 0.9871001 | 0.9882899 | 0.9930459 |
El accuracy promedio estimado usando un modelo de Regresión logística, se ha conseguido un accuracy promedio de validación del 98.71%, esto significa que el modelo predice correctamente la exclusion social por departamentos para Colombia un 98.71% de las veces.
Error de test
confusionMatrix(data = predicciones_raw_log,
reference = datos_test$hogar_exclusion,
positive = "Si")## Confusion Matrix and Statistics
##
## Reference
## Prediction Si No
## Si 2348 19
## No 40 2984
##
## Accuracy : 0.9891
## 95% CI : (0.9859, 0.9917)
## No Information Rate : 0.557
## P-Value [Acc > NIR] : < 2e-16
##
## Kappa : 0.9778
##
## Mcnemar's Test P-Value : 0.00922
##
## Sensitivity : 0.9832
## Specificity : 0.9937
## Pos Pred Value : 0.9920
## Neg Pred Value : 0.9868
## Prevalence : 0.4430
## Detection Rate : 0.4355
## Detection Prevalence : 0.4391
## Balanced Accuracy : 0.9885
##
## 'Positive' Class : Si
## # Error de test
error_test_log <- mean(predicciones_raw_log != datos_test$hogar_exclusion)
paste("El error del test:", round(error_test_log*100, 2), "%")## [1] "El error del test: 1.09 %"Modelo LDA
Ajuste, optimización y validación del modelo
Este algoritmo no tiene ningún hiperparámetro.
predicciones_lda <- predict(modelo_lda,
newdata = datos_train %>% select(-FEX_C))Una validación cruzada con 10 particiones y 5 repeticiones implica ajustar y evaluar el modelo 10 x 5 = 50 veces, cada vez con una partición distinta, y un ajuste con los datos de entrenamiento para crear el modelo final.
modelo_lda## Linear Discriminant Analysis
##
## 21566 samples
## 14 predictor
## 2 classes: 'Si', 'No'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 5 times)
## Summary of sample sizes: 19409, 19409, 19408, 19410, 19410, 19410, ...
## Resampling results:
##
## Accuracy Kappa
## 0.9805435 0.9605018modelo_lda$finalModel## Call:
## lda(x, grouping = y)
##
## Prior probabilities of groups:
## Si No
## 0.4430121 0.5569879
##
## Group means:
## control_crec icbf I_HOGAR ingresos_flia seg_alim Zona
## Si 1.183797 1.885807 -668957.5 1.347184 1.890726 1.708813
## No 1.074259 1.986930 522426.9 1.681152 1.965701 1.314186
## Registro_civil hacinnomiti quintil sgss No_Actividad_Estimulacion
## Si 1.100272 1.843521 1.721583 1.137743 1.617542
## No 1.000000 1.975191 3.986097 1.037879 1.661255
## Actividad_Deportiva Afect_Fenome_Natural DEPTO_D
## Si 1.981578 1.762822 16.62267
## No 1.954795 1.870380 16.33059
##
## Coefficients of linear discriminants:
## LD1
## control_crec -6.204096e-02
## icbf 2.022029e+00
## I_HOGAR -5.105960e-08
## ingresos_flia -9.539691e-02
## seg_alim 6.834125e-02
## Zona 2.021592e-02
## Registro_civil -1.668489e+00
## hacinnomiti 1.177529e-01
## quintil 1.285003e+00
## sgss 7.940694e-02
## No_Actividad_Estimulacion -5.786763e-02
## Actividad_Deportiva -2.467621e-01
## Afect_Fenome_Natural -7.206952e-03
## DEPTO_D -1.341810e-03summary(modelo_lda$resample$Accuracy)| Min. | 1st Qu. | Median | Mean | 3rd Qu. | Max. |
|---|---|---|---|---|---|
| 0.97405 | 0.9783239 | 0.9812196 | 0.9805435 | 0.9823829 | 0.9856282 |
El accuracy promedio estimado usando un modelo LDA, se ha conseguido un accuracy promedio de validación del 98.05%, esto significa que el modelo predice correctamente la exclusion social por departamentos para Colombia un 98.05% de las veces.
Error de test
confusionMatrix(data = predicciones_raw_lda,
reference = datos_test$hogar_exclusion,
positive = "Si")## Confusion Matrix and Statistics
##
## Reference
## Prediction Si No
## Si 2320 17
## No 68 2986
##
## Accuracy : 0.9842
## 95% CI : (0.9805, 0.9874)
## No Information Rate : 0.557
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.968
##
## Mcnemar's Test P-Value : 5.852e-08
##
## Sensitivity : 0.9715
## Specificity : 0.9943
## Pos Pred Value : 0.9927
## Neg Pred Value : 0.9777
## Prevalence : 0.4430
## Detection Rate : 0.4303
## Detection Prevalence : 0.4335
## Balanced Accuracy : 0.9829
##
## 'Positive' Class : Si
## # Error de test
error_test_lda <- mean(predicciones_raw_lda != datos_test$hogar_exclusion)
paste("El error del test:", round(error_test_lda*100, 2), "%")## [1] "El error del test: 1.58 %"Árbol de Clasificación Simple
Ajuste, optimización y validación del modelo
Este algoritmo no tiene ningún hiperparámetro.
predicciones_tree <- predict(modelo_C50Tree,
newdata = datos_train %>% select(-FEX_C))Una validación cruzada con 10 particiones y 5 repeticiones implica ajustar y evaluar el modelo 10 x 5 = 50 veces, cada vez con una partición distinta, y un ajuste con los datos de entrenamiento para crear el modelo final.
modelo_C50Tree## Single C5.0 Tree
##
## 21566 samples
## 14 predictor
## 2 classes: 'Si', 'No'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 5 times)
## Summary of sample sizes: 19409, 19409, 19408, 19410, 19410, 19410, ...
## Resampling results:
##
## Accuracy Kappa
## 0.9999073 0.9998121summary(modelo_C50Tree$finalModel)##
## Call:
## C50:::C5.0.default(x = x, y = y, weights = wts)
##
##
## C5.0 [Release 2.07 GPL Edition] Sun Jan 03 19:07:33 2021
## -------------------------------
##
## Class specified by attribute `outcome'
##
## Read 21566 cases (15 attributes) from undefined.data
##
## Decision tree:
##
## quintil <= 2: Si (8620)
## quintil > 2:
## :...Registro_civil > 1: Si (322)
## Registro_civil <= 1:
## :...icbf <= 1:
## :...No_Actividad_Estimulacion <= 1: No (91)
## : No_Actividad_Estimulacion > 1:
## : :...Actividad_Deportiva <= 1: No (66)
## : Actividad_Deportiva > 1: Si (590)
## icbf > 1:
## :...hacinnomiti > 1: No (11563/2)
## hacinnomiti <= 1:
## :...seg_alim <= 1: Si (20)
## seg_alim > 1: No (294)
##
##
## Evaluation on training data (21566 cases):
##
## Decision Tree
## ----------------
## Size Errors
##
## 8 2( 0.0%) <<
##
##
## (a) (b) <-classified as
## ---- ----
## 9552 2 (a): class Si
## 12012 (b): class No
##
##
## Attribute usage:
##
## 100.00% quintil
## 60.03% Registro_civil
## 58.54% icbf
## 55.07% hacinnomiti
## 3.46% No_Actividad_Estimulacion
## 3.04% Actividad_Deportiva
## 1.46% seg_alim
##
##
## Time: 0.1 secssummary(modelo_C50Tree$resample$Accuracy)| Min. | 1st Qu. | Median | Mean | 3rd Qu. | Max. |
|---|---|---|---|---|---|
| 0.9990724 | 1 | 1 | 0.9999073 | 1 | 1 |
El accuracy promedio estimado usando un modelo Árbol de Clasificación Simple, se ha conseguido un accuracy promedio de validación del 99.99%, esto significa que el modelo predice correctamente la exclusion social por departamentos para Colombia un 99.99% de las veces.
Error de test
confusionMatrix(data = predicciones_raw_arbol,
reference = datos_test$hogar_exclusion,
positive = "Si")## Confusion Matrix and Statistics
##
## Reference
## Prediction Si No
## Si 2387 2
## No 1 3001
##
## Accuracy : 0.9994
## 95% CI : (0.9984, 0.9999)
## No Information Rate : 0.557
## P-Value [Acc > NIR] : <2e-16
##
## Kappa : 0.9989
##
## Mcnemar's Test P-Value : 1
##
## Sensitivity : 0.9996
## Specificity : 0.9993
## Pos Pred Value : 0.9992
## Neg Pred Value : 0.9997
## Prevalence : 0.4430
## Detection Rate : 0.4428
## Detection Prevalence : 0.4431
## Balanced Accuracy : 0.9995
##
## 'Positive' Class : Si
## # Error de test
error_test_arbol <- mean(predicciones_raw_arbol != datos_test$hogar_exclusion)
paste("El error del test:", round(error_test_arbol*100, 2), "%")## [1] "El error del test: 0.06 %"Random Forest
Ajuste, optimización y validación del modelo
El método ranger de caret tiene 3 hiperparámetros:
* mtry: número predictores seleccionados aleatoriamente en cada árbol.
* min.node.size: tamaño mínimo que tiene que tener un nodo para poder ser dividido.
* splitrule: criterio de división.
predicciones_rf <- predict(modelo_rf,
newdata = datos_train %>% select(-FEX_C))Una validación cruzada con 10 particiones y 5 repeticiones implica ajustar y evaluar el modelo 10 x 5 = 50 veces, cada vez con una partición distinta, y un ajuste con los datos de entrenamiento para crear el modelo final.
modelo_rf## Random Forest
##
## 21566 samples
## 14 predictor
## 2 classes: 'Si', 'No'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 5 times)
## Summary of sample sizes: 19409, 19409, 19408, 19410, 19410, 19410, ...
## Resampling results across tuning parameters:
##
## mtry min.node.size Accuracy Kappa
## 3 2 0.9997774 0.9995489
## 3 3 0.9998053 0.9996053
## 3 4 0.9998053 0.9996053
## 3 5 0.9998145 0.9996241
## 3 10 0.9997867 0.9995677
## 3 15 0.9997125 0.9994173
## 3 20 0.9996940 0.9993797
## 3 30 0.9996012 0.9991917
## 4 2 0.9998238 0.9996429
## 4 3 0.9998145 0.9996241
## 4 4 0.9998423 0.9996805
## 4 5 0.9998331 0.9996617
## 4 10 0.9998053 0.9996053
## 4 15 0.9997682 0.9995301
## 4 20 0.9997774 0.9995489
## 4 30 0.9997311 0.9994549
## 5 2 0.9998423 0.9996805
## 5 3 0.9998423 0.9996805
## 5 4 0.9998423 0.9996805
## 5 5 0.9998331 0.9996617
## 5 10 0.9998423 0.9996805
## 5 15 0.9998331 0.9996617
## 5 20 0.9998145 0.9996241
## 5 30 0.9997589 0.9995113
## 7 2 0.9998980 0.9997933
## 7 3 0.9998887 0.9997745
## 7 4 0.9998794 0.9997557
## 7 5 0.9998980 0.9997933
## 7 10 0.9998980 0.9997933
## 7 15 0.9998702 0.9997369
## 7 20 0.9998887 0.9997745
## 7 30 0.9998887 0.9997745
##
## Tuning parameter 'splitrule' was held constant at a value of gini
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were mtry = 7, splitrule = gini
## and min.node.size = 5.modelo_rf$finalModel## Ranger result
##
## Call:
## ranger::ranger(dependent.variable.name = ".outcome", data = x, mtry = min(param$mtry, ncol(x)), min.node.size = param$min.node.size, splitrule = as.character(param$splitrule), write.forest = TRUE, probability = classProbs, ...)
##
## Type: Classification
## Number of trees: 500
## Sample size: 21566
## Number of independent variables: 14
## Mtry: 7
## Target node size: 5
## Variable importance mode: none
## Splitrule: gini
## OOB prediction error: 0.01 %ggplot(modelo_rf, highlight = TRUE) +
scale_x_continuous(breaks = 1:30) +
labs(title = "Accuracy del modelo Random Forest") +
guides(color = guide_legend(title = "mtry"),
shape = guide_legend(title = "mtry")) +
theme_bw()
Usando un modelo Random Forest con mtry = 7, min.node.size = 15 y splitrule = “gini”, se ha conseguido un accuracy promedio de validación del 99.98%.
DT::datatable(modelo_rf$resample,
extensions = list('Scroller'=TRUE ,'FixedColumns'=NULL),
options = list(autowidth=T,pageLength=20,fixedHeader=T,scrollY=200,scrollX=200,searching = T,
fixedColumns=list(leftColumns=1) )) %>% formatPercentage(1:2, 3)summary(modelo_rf$resample$Accuracy)| Min. | 1st Qu. | Median | Mean | 3rd Qu. | Max. |
|---|---|---|---|---|---|
| 0.9972171 | 0.9995364 | 1 | 0.999816 | 1 | 1 |
El accuracy promedio estimado usando un modelo Random Forest, se ha conseguido un accuracy promedio de validación del 99.98%, esto significa que el modelo predice correctamente la exclusion social por departamentos para Colombia un 99.98% de las veces.
Error de test
confusionMatrix(data = predicciones_raw_rf,
reference = datos_test$hogar_exclusion,
positive = "Si")## Confusion Matrix and Statistics
##
## Reference
## Prediction Si No
## Si 2387 2
## No 1 3001
##
## Accuracy : 0.9994
## 95% CI : (0.9984, 0.9999)
## No Information Rate : 0.557
## P-Value [Acc > NIR] : <2e-16
##
## Kappa : 0.9989
##
## Mcnemar's Test P-Value : 1
##
## Sensitivity : 0.9996
## Specificity : 0.9993
## Pos Pred Value : 0.9992
## Neg Pred Value : 0.9997
## Prevalence : 0.4430
## Detection Rate : 0.4428
## Detection Prevalence : 0.4431
## Balanced Accuracy : 0.9995
##
## 'Positive' Class : Si
## # Error de test
error_test_rf <- mean(predicciones_raw_rf != datos_test$hogar_exclusion)
paste("El error del test:", round(error_test_rf*100, 2), "%")## [1] "El error del test: 0.06 %"Comparando los modelos
Una vez que se han entrenado y optimizado distintos modelos, se tiene que identificar cuál de ellos consigue mejores resultados para el problema en cuestión, en este caso, predecir la presencia de exclusion social en los hogares por cada uno de los departamentos Colombianos.
Métricas de validación
modelos <- list(KNN = modelo_knn, logistic = modelo_logistic,
lda= modelo_lda, arbol = modelo_C50Tree, rf = modelo_rf)
resultados_resamples <- resamples(modelos)
DT::datatable(resultados_resamples$values,
extensions = list('Scroller'=TRUE ,'FixedColumns'=NULL),
options = list(autowidth=T,pageLength=20,fixedHeader=T,scrollY=200,scrollX=200,searching = T,
fixedColumns=list(leftColumns=1) ))%>% formatPercentage(2:11, 3)metricas_resamples <- resultados_resamples$values %>%
gather(key = "modelo", value = "valor", -Resample) %>%
separate(col = "modelo", into = c("modelo", "metrica"),
sep = "~", remove = TRUE)DT::datatable(metricas_resamples,
extensions = list('Scroller'=TRUE ,'FixedColumns'=NULL),
options = list(autowidth=T,pageLength=20,fixedHeader=T,scrollY=200,scrollX=200,searching = T,
fixedColumns=list(leftColumns=1) ))%>% formatPercentage('valor', 3)plot_metricas## Warning: Removed 43 rows containing missing values (geom_point).## Warning: Removed 1 rows containing missing values (geom_text).
Para determinar si las diferencias entre los modelos generados son significativas, se recurre a test estadísticos.
Test de Friedman para comparar el accuracy de los modelos
matriz_metricas <- metricas_resamples %>% filter(metrica == "Accuracy") %>%
spread(key = modelo, value = valor) %>%
select(-Resample, -metrica) %>% as.matrix()
friedman.test(y = matriz_metricas)##
## Friedman rank sum test
##
## data: matriz_metricas
## Friedman chi-squared = 199.79, df = 4, p-value < 2.2e-16Para un nivel de significancia (α = 0.05), el test de Friedman sí encuentra evidencias para rechazar la hipótesis nula de que los 5 clasificadores consiguen la misma precisión.
kable(comparaciones, digits = 3, row.names = FALSE, align = "c",caption = NULL)| modeloA | modeloB | p_value |
|---|---|---|
| KNN | arbol | 0 |
| lda | arbol | 0 |
| lda | KNN | 0 |
| logistic | arbol | 0 |
| logistic | KNN | 0 |
| logistic | lda | 0 |
| rf | arbol | 1 |
| rf | KNN | 0 |
| rf | lda | 0 |
| rf | logistic | 0 |
Acorde a las comparaciones por pares, no existen evidencias suficientes para considerar que la capacidad predictiva de los modelos Random Forest y Arbol de Clasificación es distinta.
Error de test
A continuación se muestran las predicciones de cada uno de los modelos, tanto las observaciones de entrenamiento como para las de test. Con toda esta información, se compararan los resultados de predicción entre modelos y las diferencias entre conjunto de entrenamiento y test.
datos_test <- select(datos_test, -Sexo, -cuidado_menor,-Estado_de_Salud,-Victima_Hecho_Viol,-vacunas,-Edad)
testY <- datos_test %>% select(-FEX_C)
predicciones <- extractPrediction(
models = modelos,
testX = datos_test %>% select(-hogar_exclusion,-FEX_C),
testY = testY$hogar_exclusion )DT::datatable(predicciones,
extensions = list('Scroller'=TRUE ,'FixedColumns'=NULL),
options = list(autowidth=T,pageLength=20,fixedHeader=T,scrollY=200,scrollX=200,searching = T,
fixedColumns=list(leftColumns=1) ))## Warning in instance$preRenderHook(instance): It seems your data is too big
## for client-side DataTables. You may consider server-side processing: https://
## rstudio.github.io/DT/server.htmlCon toda esta información, es sencillo comparar los resultados de predicción entre modelos y las diferencias entre conjunto de entrenamiento y test.
metricas_predicciones$accuracy## [1] 0.9850799Puede verse que, todos los modelos, consiguen más predicciones correctas en el conjunto de entrenamiento que en el de test, de ahí que las métricas obtenidas en el entrenamiento no deban utilizarse para evaluar los modelos, son excesivamente optimistas. El modelo arbol y rf consiguen el accuracy de test más alto.
Cuando he comparado los modelos podria concluir que el mejor modelo es el basado en Arboles Simples y el Random Forest.
comparaciones %>% filter((modeloA == "rf") | (modeloA == "arbol" & modeloB == "rf"))| modeloA | modeloB | p_value |
|---|---|---|
| rf | arbol | 1 |
| rf | KNN | 0 |
| rf | lda | 0 |
| rf | logistic | 0 |
7. Análisis Clúster
Análisis Clúster por Departamentos
Se ha usado el método del codo o elbow method para estimar el número K óptimo de clusters cuando no se dispone de información adicional en la que basarse, es aplicar el algoritmo de K-means para un rango de valores de K e identificar aquel valor a partir del cual la reducción en la suma total de varianza intra-cluster deja de ser sustancial.
Ahora continuamos nuestro Analisis con los Datos Test, y confirmar que han sido bien clasificados deacuerdo con la realidad del pais (seleccionamos las predicciones del arbol simple).
fviz_nbclust(x = survival, FUNcluster = kmeans, method = "wss", k.max = 15,
diss = get_dist(survival, method = "euclidean"), nstart = 50)
En este caso el K óptimo es 3.
set.seed(123)
km_clusters <- kmeans(x = survival,
centers = 3,
nstart = 50)
fviz_cluster(object = km_clusters, data = survival, show.clust.cent = TRUE,
ellipse.type = "euclid", star.plot = TRUE, repel = TRUE) +
labs(title = "Resultados clustering K-means") +
theme_bw() +
theme(legend.position = "none")
En el gráfico anterior se puede visualizar las agrupaciones resultantes del clúster. Si el número de variables (dimensionalidad) es mayor de 2, automáticamente realiza un PCA y representa las dos primeras componentes principales.
Al aplicar hierarchical clustering, empleando como medida de similitud la distancia euclídea y linkage complete, se obtiene el siguiente dendrograma.
set.seed(101)
hc_euclidea_completo <- hclust(d = dist(x = survival, method = "euclidean"),
method = "complete")
fviz_dend(x = hc_euclidea_completo, k = 3, cex = 0.4) +
geom_hline(yintercept = 0.65, linetype = "dashed") +
labs(title = "Herarchical clustering",
subtitle = "Distancia euclídeana, Lincage complete, K=3")
En la base del dendrograma, cada observación forma una terminación individual conocida como hoja o leaf del árbol (clasificación de los departamentos).
Análsis Clúster por Variables
# Transpone todas las columnas menos la primer
df_transpose <- data.frame(t(survival[]))
df_transpose <- df_transpose[ !(rownames(df_transpose) %in% c("hogar_exclusionSi", "predicciones_raw_arbolSi")), ]Al aplicar hierarchical clustering, empleando como medida de similitud la distancia euclídea y linkage complete, se obtiene el siguiente dendrograma para las variables.
set.seed(101)
hc_euclidea_completo <- hclust(d = dist(x = df_transpose, method = "euclidean"),
method = "complete")
fviz_dend(x = hc_euclidea_completo, k = 3, cex = 0.4) +
geom_hline(yintercept = 1.5, linetype = "dashed") +
labs(title = "Herarchical clustering",
subtitle = "Distancia euclídeana, Lincage complete, K=3")
Una forma menos frecuente de representar los resultados de un hierarchical clustering es combinándolos con una reducción de dimensionalidad por PCA. Primero, se calculan las componentes principales y se representan las observaciones en un scatterplot empleando las dos primeras componentes, finalmente se colorean los clusters mediante elipses.
fviz_cluster(object = list(data=df_transpose, cluster=cutree(hc_euclidea_completo, k=3)),
ellipse.type = "convex", repel = TRUE, show.clust.cent = FALSE,
labelsize = 8) +
labs(title = "Hierarchical clustering + Proyección PCA",
subtitle = "Distancia euclídeana, Lincage complete, K=3") +
theme_bw() +
theme(legend.position = "bottom")
8. Indicador de Exclusión Social
De acuerdo con lo explicado anteriormente se presenta a continuación un resumen del escenario en donde se ve por departamentos cual es el porcentaje de hogares encuestados que cuentan o están en riesgo de exclusión social de acuerdo con las características del hogar encuestado.
Se hace un resumen a su vez por cada una de las variables seleccionadas de acuerdo con sus características.
indicadorMapa Geopolitico A continuación se muestra uno de los mapa geopoliticos Colombianos creados en el software ArcMap 10.5 (ArcGIS), en el que se puede ver la distribución de la exclusión social por departamentos para Colombia y que haciendo una mirada critica se asemeja claramente a la realidad del país. (El resto de los mapas se podran ver en los anexos de la misma).
knitr::include_graphics('hogar_exclusion.PNG')
Información sesión
sesion_info <- devtools::session_info()
librerias <- dplyr::select(
tibble::as_tibble(sesion_info$packages),
c(package, loadedversion, source) )
head(librerias)| package | loadedversion | source |
|---|---|---|
| abind | 1.4-5 | CRAN (R 4.0.0) |
| assertthat | 0.2.1 | CRAN (R 4.0.2) |
| backports | 1.1.9 | CRAN (R 4.0.2) |
| base64enc | 0.1-3 | CRAN (R 4.0.0) |
| bitops | 1.0-6 | CRAN (R 4.0.0) |
| blob | 1.2.1 | CRAN (R 4.0.2) |